Repeated Rarefaction Function

Load in neccesary data

#sample_info <- read.table("sample_data.txt")
sample_info_tab <- read.table("/Users/magdalenapolder/Documents/examensarbete/scripting/sample_info.tsv", header=T, row.names=1, check.names=F, sep="\t")[-28, ]
sample_info_tab$color <- as.character(sample_info_tab$color)
count_data <- read.table("ASVs_counts.tsv", header=T, row.names=1,
             check.names=F, sep="\t")

Load neccesary packages

library(fossil)

Loading required package: sp
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, were retired in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.

Attaching package: ‘sp’

The following object is masked from ‘package:IRanges’:

    %over%

Loading required package: maps

Attaching package: ‘maps’

The following object is masked from ‘package:purrr’:

    map

Loading required package: shapefiles
Loading required package: foreign

Attaching package: ‘shapefiles’

The following objects are masked from ‘package:foreign’:

    read.dbf, write.dbf

Function

#Define Input
count_input <- count_data
info_input <- sample_info_tab
rarefaction_threshhold <- 50
repeat_amount <- 1


#Set up working files
rarified_count <- rrarefy(t(count_input),rarefaction_threshhold)
duplicated_info <- info_input

#Perform repeated rarefaction
if (repeat_amount > 1) {
  for (x in 2:repeat_amount){
    rarified_count <- rbind(rarified_count,rrarefy(t(count_input),2000))
    duplicated_info <- rbind(duplicated_info, info_input)
    }
  }



rare_count_phy_repeat <- otu_table(t(rarified_count), taxa_are_rows=T)

sample_info_tab_phy_repeat <- sample_data(duplicated_info)
rare_physeq_repeat <- phyloseq(rare_count_phy_repeat,sample_info_tab_phy_repeat)

vst_pcoa_repeat <- ordinate(rare_physeq_repeat, method="PCoA", distance="bray")

plot_ordination(rare_physeq_repeat, vst_pcoa_repeat,  color = "location")

coordinates <- vst_pcoa_repeat$vectors[,c(1,2)]
sum_test <- as.data.frame(rowsum(coordinates, row.names(coordinates)))
scale_factor <- length(coordinates)/length(t(sum_test))
sum_test <- as.data.frame(sum_test/scale_factor)

my_plot <- (plot_ordination(rare_physeq_repeat, vst_pcoa_repeat,justDF = TRUE))

mean_df <- my_plot[1:length(sum_test$Axis.1),]
mean_df[,1] <- sum_test$Axis.1
mean_df[,2] <- sum_test$Axis.2

#plot(x=sum_test[,1], y=sum_test[,2])

ggplot(data = my_plot, aes(x=my_plot[,1], y=my_plot[,2],color=my_plot$location)) + labs(colour = "Location", x = "NMDS1", y ="NMDS2") + geom_point(na.rm=TRUE, shape=NA) +stat_ellipse(linetype = 1,lwd = 0.8, aes(color=my_plot$location, group=my_plot$sample_id))+geom_point(data=mean_df, mapping =aes(x=Axis.1, y=Axis.2, alpha=0, color=mean_df$location)) + geom_point(alpha = 0, na.rm=TRUE)

#plot_ordination(rare_physeq_repeat, vst_pcoa_repeat,) + geom_point(shape="cross")+
#  geom_point(data=sum_test, mapping =aes(x=MDS1, y=MDS2)) +
#  stat_ellipse(linetype = 1,lwd = 0.8) #+
  #geom_point(size=1) + labs(col="location") + 
    #geom_text(aes(label=rownames(duplicated_info), hjust=0.3, vjust=-0.4), size=3)

repeat_raref <- function(count, info, repeats, threshold, method, colorb, shapeb, cloud = FALSE, ellipse = TRUE) {
  #Set up working files
  rarified_count <- rrarefy(t(count),threshold)
  duplicated_info <- info
  
  #Perform repeated rarefaction
  if (repeats > 1) {
    for (x in 2:repeats){
      rarified_count <- rbind(rarified_count,rrarefy(t(count_input),2000))
      duplicated_info <- rbind(duplicated_info, info)
      }
    }

  #Convert the input into physeq objects
  rare_count_phy_repeat <- otu_table(t(rarified_count), taxa_are_rows=T)
  sample_info_tab_phy_repeat <- sample_data(duplicated_info)
  rare_physeq_repeat <- phyloseq(rare_count_phy_repeat,sample_info_tab_phy_repeat)
  
  #Perform the Ordination calculation
  vst_pcoa <- ordinate(rare_physeq_repeat, method=method, distance="bray")
  
  if (method == "NMDS") {
    xaxis <- "NMDS1"
    yaxis <- "NMDS2"
    #Create matrix with mean position for each sample
    coordinates <- vst_pcoa$points
    sum_test <- as.data.frame(rowsum(coordinates, row.names(coordinates)))
    scale_factor <- length(coordinates)/length(t(sum_test))
    sum_test <- as.data.frame(sum_test/scale_factor)
  
    #Convert ordination result into a data frame object
    my_plot <- (plot_ordination(rare_physeq_repeat, vst_pcoa,justDF = TRUE))
  
    #Add info to mean location data frame
    mean_df <- my_plot[1:length(sum_test$MDS1),]
    mean_df[,1] <- sum_test$MDS1
    mean_df[,2] <- sum_test$MDS2
  }
  
  if (method == "PCoA") {
    xaxis <- "Axis.1"
    yaxis <- "Axis.2"
    #Create matrix with mean position for each sample
    coordinates <- vst_pcoa$vectors[,c(1,2)]
    sum_test <- as.data.frame(rowsum(coordinates, row.names(coordinates)))
    scale_factor <- length(coordinates)/length(t(sum_test))
    sum_test <- as.data.frame(sum_test/scale_factor)
    
    my_plot <- (plot_ordination(rare_physeq_repeat, vst_pcoa,justDF = TRUE))
    
    mean_df <- my_plot[1:length(sum_test$Axis.1),]
    mean_df[,1] <- sum_test$Axis.1
    mean_df[,2] <- sum_test$Axis.2
  }
  
  
  
  color_mean <- unlist(mean_df[colorb])
  color_tot <- unlist(my_plot[colorb])
  
  if (!is.na(shapeb)) {
    shape_mean <- unlist(mean_df[shapeb])
    shape_tot <- unlist(my_plot[shapeb])
  }
  
  if (repeats <= 3 && ellipse==TRUE){
    print("Not printing ellipse due to too few repeats. Need at least 4 repeats to calculate confidence intervals.")
    ellipse = FALSE
  }
  
  #Create the plot and print it
  finished_plot <- ggplot(data = my_plot, aes(x=my_plot[,1], y=my_plot[,2],color=color_tot, shape=shape_tot)) + 
    labs(colour = colorb, shape = shapeb, x = xaxis, y =yaxis) + 
    {if(cloud)geom_point(na.rm=TRUE)} +
    {if(ellipse)stat_ellipse(linetype = 1,lwd = 0.8, aes(color=color_tot, group=my_plot$sample_id))} +
    {if(!cloud)geom_point(data=mean_df, mapping =aes(x=mean_df[,1], y=mean_df[,2], color=color_mean, shape=shape_mean))}
  
  print(finished_plot)
  return(my_plot)
}

complete_ordination <- repeat_raref(count_data, sample_info_tab, repeats = 10, threshold = 1000, method = "NMDS","location", "type", cloud= TRUE, ellipse=FALSE)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.2086504 
Run 1 stress 0.2563442 
Run 2 stress 0.2259494 
Run 3 stress 0.2575658 
Run 4 stress 0.2283044 
Run 5 stress 0.245765 
Run 6 stress 0.2486569 
Run 7 stress 0.245768 
Run 8 stress 0.2378086 
Run 9 stress 0.2490392 
Run 10 stress 0.2516118 
Run 11 stress 0.2231128 
Run 12 stress 0.2612058 
Run 13 stress 0.2116076 
Run 14 stress 0.2266804 
Run 15 stress 0.2490038 
Run 16 stress 0.2093084 
Run 17 stress 0.2564273 
Run 18 stress 0.2093084 
Run 19 stress 0.2092393 
Run 20 stress 0.2396593 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
     6: stress ratio > sratmax
    11: scale factor of the gradient < sfgrmin

repeat_list <- c(10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 80, 90, 100, 125, 150, 175, 200)
time_data <- matrix(ncol = 4, nrow = 0)
colnames(time_data) <- c("Repeat Amount", "Threshold", "Time in sec", "Time in min" )


for (x in repeat_list) {
  print(paste("Running with ",x," repeats"))
  time1 <- Sys.time()
  complete_ordination <- repeat_raref(count_data, sample_info_tab, repeats = x, threshold = 1000, "sample_id", "location")
  time2 <-Sys.time()
  time_taken_secs <- difftime(time2, time1, units="secs")
  time_taken_min <- difftime(time2, time1, units="mins")
  time_data <- rbind(time_data, c(x, 1000, time_taken_secs, time_taken_min))
}

[1] "Running with  10  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2068946 
Run 1 stress 0.2663137 
Run 2 stress 0.208834 
Run 3 stress 0.2108935 
Run 4 stress 0.2484204 
Run 5 stress 0.2472172 
Run 6 stress 0.2341156 
Run 7 stress 0.2587074 
Run 8 stress 0.2148219 
Run 9 stress 0.2122492 
Run 10 stress 0.2671656 
Run 11 stress 0.2366999 
Run 12 stress 0.2394257 
Run 13 stress 0.2363499 
Run 14 stress 0.2392263 
Run 15 stress 0.2484121 
Run 16 stress 0.2092305 
Run 17 stress 0.2329405 
Run 18 stress 0.237903 
Run 19 stress 0.2287428 
Run 20 stress 0.2107879 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
    12: stress ratio > sratmax
     5: scale factor of the gradient < sfgrmin
[1] "Running with  15  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2037621 
Run 1 stress 0.2144868 
Run 2 stress 0.2339792 
Run 3 stress 0.2076302 
Run 4 stress 0.266159 
Run 5 stress 0.2037621 
... Procrustes: rmse 1.353935e-05  max resid 8.49778e-05 
... Similar to previous best
Run 6 stress 0.2620269 
Run 7 stress 0.2037621 
... Procrustes: rmse 1.964777e-05  max resid 0.0001768671 
... Similar to previous best
Run 8 stress 0.2278356 
Run 9 stress 0.2531525 
Run 10 stress 0.2219272 
Run 11 stress 0.2404791 
Run 12 stress 0.2309629 
Run 13 stress 0.2050069 
Run 14 stress 0.2134233 
Run 15 stress 0.2485466 
Run 16 stress 0.2721227 
Run 17 stress 0.2063026 
Run 18 stress 0.2062653 
Run 19 stress 0.2663672 
Run 20 stress 0.2544041 
*** Best solution repeated 2 times
[1] "Running with  20  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2172838 
Run 1 stress 0.247867 
Run 2 stress 0.2323725 
Run 3 stress 0.215454 
... New best solution
... Procrustes: rmse 0.008543489  max resid 0.07234141 
Run 4 stress 0.2104664 
... New best solution
... Procrustes: rmse 0.009625521  max resid 0.06633488 
Run 5 stress 0.2324129 
Run 6 stress 0.2582749 
Run 7 stress 0.2092665 
... New best solution
... Procrustes: rmse 0.00585217  max resid 0.07313217 
Run 8 stress 0.2102615 
Run 9 stress 0.2482576 
Run 10 stress 0.2107992 
Run 11 stress 0.2062045 
... New best solution
... Procrustes: rmse 0.005213015  max resid 0.05727515 
Run 12 stress 0.24784 
Run 13 stress 0.2556115 
Run 14 stress 0.2319646 
Run 15 stress 0.2621871 
Run 16 stress 0.2278621 
Run 17 stress 0.2113299 
Run 18 stress 0.2736252 
Run 19 stress 0.2080111 
Run 20 stress 0.2081056 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
     4: stress ratio > sratmax
    13: scale factor of the gradient < sfgrmin
[1] "Running with  25  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2159204 
Run 1 stress 0.2109334 
... New best solution
... Procrustes: rmse 0.01345458  max resid 0.0647212 
Run 2 stress 0.2323745 
Run 3 stress 0.2096228 
... New best solution
... Procrustes: rmse 0.00694352  max resid 0.05694282 
Run 4 stress 0.2302605 
Run 5 stress 0.2333996 
Run 6 stress 0.2242031 
Run 7 stress 0.2061253 
... New best solution
... Procrustes: rmse 0.00675013  max resid 0.07979637 
Run 8 stress 0.2481123 
Run 9 stress 0.2386248 
Run 10 stress 0.2131409 
Run 11 stress 0.2100337 
Run 12 stress 0.2562804 
Run 13 stress 0.2039217 
... New best solution
... Procrustes: rmse 0.004461687  max resid 0.03551654 
Run 14 stress 0.2133143 
Run 15 stress 0.2050656 
Run 16 stress 0.2322626 
Run 17 stress 0.2455612 
Run 18 stress 0.2270098 
Run 19 stress 0.2324001 
Run 20 stress 0.2292402 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: scale factor of the gradient < sfgrmin
[1] "Running with  30  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2048868 
Run 1 stress 0.2066875 
Run 2 stress 0.2502874 
Run 3 stress 0.2077659 
Run 4 stress 0.2086475 
Run 5 stress 0.227905 
Run 6 stress 0.2256868 
Run 7 stress 0.2322326 
Run 8 stress 0.2456374 
Run 9 stress 0.2153393 
Run 10 stress 0.235093 
Run 11 stress 0.2416219 
Run 12 stress 0.2322977 
Run 13 stress 0.2491814 
Run 14 stress 0.2123347 
Run 15 stress 0.2081935 
Run 16 stress 0.2073044 
Run 17 stress 0.2079442 
Run 18 stress 0.2427943 
Run 19 stress 0.2475072 
Run 20 stress 0.2485183 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
     2: stress ratio > sratmax
    15: scale factor of the gradient < sfgrmin
[1] "Running with  35  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2159627 
Run 1 stress 0.2440572 
Run 2 stress 0.2325243 
Run 3 stress 0.2165349 
Run 4 stress 0.2365486 
Run 5 stress 0.2434068 
Run 6 stress 0.2120876 
... New best solution
... Procrustes: rmse 0.01209213  max resid 0.05763703 
Run 7 stress 0.2403782 
Run 8 stress 0.2415735 
Run 9 stress 0.2331229 
Run 10 stress 0.2611486 
Run 11 stress 0.2381923 
Run 12 stress 0.215076 
Run 13 stress 0.2491559 
Run 14 stress 0.2372107 
Run 15 stress 0.2097071 
... New best solution
... Procrustes: rmse 0.005385369  max resid 0.0493466 
Run 16 stress 0.2184571 
Run 17 stress 0.2646963 
Run 18 stress 0.2270707 
Run 19 stress 0.2363402 
Run 20 stress 0.2656087 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
     1: stress ratio > sratmax
    14: scale factor of the gradient < sfgrmin
[1] "Running with  40  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2148911 
Run 1 stress 0.2341205 
Run 2 stress 0.260825 
Run 3 stress 0.2493988 
Run 4 stress 0.2370208 
Run 5 stress 0.2376256 
Run 6 stress 0.2613994 
Run 7 stress 0.2419957 
Run 8 stress 0.2400575 
Run 9 stress 0.2500377 
Run 10 stress 0.2500017 
Run 11 stress 0.248826 
Run 12 stress 0.2099825 
... New best solution
... Procrustes: rmse 0.01046906  max resid 0.05603305 
Run 13 stress 0.2580906 
Run 14 stress 0.2055665 
... New best solution
... Procrustes: rmse 0.004565162  max resid 0.04443016 
Run 15 stress 0.246378 
Run 16 stress 0.2048284 
... New best solution
... Procrustes: rmse 0.003135009  max resid 0.04922089 
Run 17 stress 0.2309908 
Run 18 stress 0.2343338 
Run 19 stress 0.2085154 
Run 20 stress 0.2419515 
*** Best solution was not repeated -- monoMDS stopping criteria:
     7: no. of iterations >= maxit
    13: scale factor of the gradient < sfgrmin
[1] "Running with  45  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2155541 
Run 1 stress 0.2456793 
Run 2 stress 0.2422816 
Run 3 stress 0.2253852 
Run 4 stress 0.2374882 
Run 5 stress 0.2448256 
Run 6 stress 0.225317 
Run 7 stress 0.25378 
Run 8 stress 0.2476409 
Run 9 stress 0.2128044 
... New best solution
... Procrustes: rmse 0.01078795  max resid 0.05583074 
Run 10 stress 0.2105192 
... New best solution
... Procrustes: rmse 0.005344977  max resid 0.04347715 
Run 11 stress 0.2049319 
... New best solution
... Procrustes: rmse 0.004572584  max resid 0.06137488 
Run 12 stress 0.2112032 
Run 13 stress 0.2445934 
Run 14 stress 0.2138173 
Run 15 stress 0.2631151 
Run 16 stress 0.2447748 
Run 17 stress 0.2468428 
Run 18 stress 0.2385953 
Run 19 stress 0.2643583 
Run 20 stress 0.205687 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
    15: scale factor of the gradient < sfgrmin
[1] "Running with  50  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2146711 
Run 1 stress 0.247094 
Run 2 stress 0.2468967 
Run 3 stress 0.225205 
Run 4 stress 0.2373683 
Run 5 stress 0.206678 
... New best solution
... Procrustes: rmse 0.008840808  max resid 0.05044393 
Run 6 stress 0.2625275 
Run 7 stress 0.2646484 
Run 8 stress 0.2435206 
Run 9 stress 0.2123092 
Run 10 stress 0.2588298 
Run 11 stress 0.2076638 
Run 12 stress 0.2057752 
... New best solution
... Procrustes: rmse 0.003290995  max resid 0.044621 
Run 13 stress 0.2624739 
Run 14 stress 0.2124649 
Run 15 stress 0.2041138 
... New best solution
... Procrustes: rmse 0.002811395  max resid 0.04392208 
Run 16 stress 0.227348 
Run 17 stress 0.2149618 
Run 18 stress 0.2481491 
Run 19 stress 0.2539749 
Run 20 stress 0.2066677 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: scale factor of the gradient < sfgrmin
[1] "Running with  55  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2149181 
Run 1 stress 0.2094504 
... New best solution
... Procrustes: rmse 0.009131158  max resid 0.04576863 
Run 2 stress 0.2630586 
Run 3 stress 0.2458759 
Run 4 stress 0.2291794 
Run 5 stress 0.2083754 
... New best solution
... Procrustes: rmse 0.004126547  max resid 0.03881932 
Run 6 stress 0.231999 
Run 7 stress 0.2452261 
Run 8 stress 0.2456985 
Run 9 stress 0.2229702 
Run 10 stress 0.2455323 
Run 11 stress 0.2602408 
Run 12 stress 0.2468942 
Run 13 stress 0.4208333 
Run 14 stress 0.2473179 
Run 15 stress 0.2249847 
Run 16 stress 0.2313205 
Run 17 stress 0.2249866 
Run 18 stress 0.2061198 
... New best solution
... Procrustes: rmse 0.004316067  max resid 0.03895677 
Run 19 stress 0.2457703 
Run 20 stress 0.2578304 
*** Best solution was not repeated -- monoMDS stopping criteria:
    20: scale factor of the gradient < sfgrmin
[1] "Running with  60  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2156977 
Run 1 stress 0.2096811 
... New best solution
... Procrustes: rmse 0.007936442  max resid 0.04992578 
Run 2 stress 0.2115546 
Run 3 stress 0.4208715 
Run 4 stress 0.2078592 
... New best solution
... Procrustes: rmse 0.003649462  max resid 0.03922537 
Run 5 stress 0.2534319 
Run 6 stress 0.2261133 
Run 7 stress 0.2467298 
Run 8 stress 0.2596728 
Run 9 stress 0.2230132 
Run 10 stress 0.4208726 
Run 11 stress 0.2081206 
... Procrustes: rmse 0.00344311  max resid 0.03870427 
Run 12 stress 0.2118499 
Run 13 stress 0.2595362 
Run 14 stress 0.2571937 
Run 15 stress 0.2286469 
Run 16 stress 0.2403542 
Run 17 stress 0.2302645 
Run 18 stress 0.2134467 
Run 19 stress 0.2645 
Run 20 stress 0.2048229 
... New best solution
... Procrustes: rmse 0.002946028  max resid 0.03562683 
*** Best solution was not repeated -- monoMDS stopping criteria:
     8: no. of iterations >= maxit
    12: scale factor of the gradient < sfgrmin
[1] "Running with  65  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.215196 
Run 1 stress 0.2367409 
Run 2 stress 0.2367428 
Run 3 stress 0.2102159 
... New best solution
... Procrustes: rmse 0.008063984  max resid 0.04271201 
Run 4 stress 0.2454367 
Run 5 stress 0.205272 
... New best solution
... Procrustes: rmse 0.003309468  max resid 0.0381221 
Run 6 stress 0.2078369 
Run 7 stress 0.2299384 
Run 8 stress 0.2068666 
Run 9 stress 0.2042599 
... New best solution
... Procrustes: rmse 0.001529567  max resid 0.03292068 
Run 10 stress 0.2329084 
Run 11 stress 0.21946 
Run 12 stress 0.2469003 
Run 13 stress 0.2446621 
Run 14 stress 0.2043244 
... Procrustes: rmse 0.0004486976  max resid 0.01891409 
Run 15 stress 0.2090608 
Run 16 stress 0.2450009 
Run 17 stress 0.2489648 
Run 18 stress 0.2289682 
Run 19 stress 0.2244397 
Run 20 stress 0.2480431 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
    15: scale factor of the gradient < sfgrmin
[1] "Running with  70  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2156611 
Run 1 stress 0.2463472 
Run 2 stress 0.2547421 
Run 3 stress 0.4209413 
Run 4 stress 0.2236026 
Run 5 stress 0.2281274 
Run 6 stress 0.2575193 
Run 7 stress 0.2308123 
Run 8 stress 0.2100399 
... New best solution
... Procrustes: rmse 0.007666282  max resid 0.04136695 
Run 9 stress 0.2082502 
... New best solution
... Procrustes: rmse 0.003740293  max resid 0.03859532 
Run 10 stress 0.2437089 
Run 11 stress 0.2635594 
Run 12 stress 0.2112224 
Run 13 stress 0.2447727 
Run 14 stress 0.2377362 
Run 15 stress 0.2593871 
Run 16 stress 0.2101711 
Run 17 stress 0.2107413 
Run 18 stress 0.2121219 
Run 19 stress 0.242272 
Run 20 stress 0.2107429 
*** Best solution was not repeated -- monoMDS stopping criteria:
     4: no. of iterations >= maxit
    16: scale factor of the gradient < sfgrmin
[1] "Running with  80  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.215822 
Run 1 stress 0.2071136 
... New best solution
... Procrustes: rmse 0.006930798  max resid 0.03995484 
Run 2 stress 0.2069358 
... New best solution
... Procrustes: rmse 0.002463907  max resid 0.03575849 
Run 3 stress 0.2463148 
Run 4 stress 0.2082768 
Run 5 stress 0.2606677 
Run 6 stress 0.2110424 
Run 7 stress 0.2284578 
Run 8 stress 0.2593379 
Run 9 stress 0.258091 
Run 10 stress 0.2426367 
Run 11 stress 0.2084672 
Run 12 stress 0.263565 
Run 13 stress 0.2610747 
Run 14 stress 0.2444638 
Run 15 stress 0.2597474 
Run 16 stress 0.2257589 
Run 17 stress 0.2482234 
Run 18 stress 0.2488592 
Run 19 stress 0.2470879 
Run 20 stress 0.2466732 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: scale factor of the gradient < sfgrmin
[1] "Running with  90  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2035858 
Run 1 stress 0.2376481 
Run 2 stress 0.2328945 
Run 3 stress 0.2116773 
Run 4 stress 0.2071597 
Run 5 stress 0.2318099 
Run 6 stress 0.2320072 
Run 7 stress 0.421032 
Run 8 stress 0.2457167 
Run 9 stress 0.2463062 
Run 10 stress 0.2441054 
Run 11 stress 0.2226625 
Run 12 stress 0.421033 
Run 13 stress 0.2628761 
Run 14 stress 0.2469244 
Run 15 stress 0.2448891 
Run 16 stress 0.2503415 
Run 17 stress 0.2460393 
Run 18 stress 0.421027 
Run 19 stress 0.245976 
Run 20 stress 0.2427347 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: scale factor of the gradient < sfgrmin
[1] "Running with  100  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2044245 
Run 1 stress 0.2370089 
Run 2 stress 0.2327077 
Run 3 stress 0.2119449 
Run 4 stress 0.2439213 
Run 5 stress 0.2341286 
Run 6 stress 0.2093629 
Run 7 stress 0.208026 
Run 8 stress 0.2132467 
Run 9 stress 0.2102858 
Run 10 stress 0.2644138 
Run 11 stress 0.2141068 
Run 12 stress 0.2472585 
Run 13 stress 0.240561 
Run 14 stress 0.2248437 
Run 15 stress 0.2332474 
Run 16 stress 0.2604631 
Run 17 stress 0.2187768 
Run 18 stress 0.2247018 
Run 19 stress 0.2333012 
Run 20 stress 0.2134549 
*** Best solution was not repeated -- monoMDS stopping criteria:
     7: no. of iterations >= maxit
    13: scale factor of the gradient < sfgrmin
[1] "Running with  125  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2159332 
Run 1 stress 0.246082 
Run 2 stress 0.2480258 
Run 3 stress 0.2470711 
Run 4 stress 0.2294322 
Run 5 stress 0.2391229 
Run 6 stress 0.23303 
Run 7 stress 0.204853 
... New best solution
... Procrustes: rmse 0.005227148  max resid 0.03211169 
Run 8 stress 0.4211179 
Run 9 stress 0.2302452 
Run 10 stress 0.2444754 
Run 11 stress 0.2451402 
Run 12 stress 0.2119602 
Run 13 stress 0.2474627 
Run 14 stress 0.2281534 
Run 15 stress 0.2475207 
Run 16 stress 0.232335 
Run 17 stress 0.2647527 
Run 18 stress 0.2402844 
Run 19 stress 0.2322469 
Run 20 stress 0.2285175 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
    17: scale factor of the gradient < sfgrmin
[1] "Running with  150  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2157534 
Run 1 stress 0.246089 
Run 2 stress 0.2475034 
Run 3 stress 0.4211591 
Run 4 stress 0.2330081 
Run 5 stress 0.2360439 
Run 6 stress 0.2357407 
Run 7 stress 0.2268399 
Run 8 stress 0.2105982 
... New best solution
... Procrustes: rmse 0.005350735  max resid 0.03034189 
Run 9 stress 0.2587474 
Run 10 stress 0.2255851 
Run 11 stress 0.2641206 
Run 12 stress 0.2442011 
Run 13 stress 0.226946 
Run 14 stress 0.2382041 
Run 15 stress 0.4211571 
Run 16 stress 0.2261881 
Run 17 stress 0.2331532 
Run 18 stress 0.4211552 
Run 19 stress 0.2447606 
Run 20 stress 0.263867 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
    15: scale factor of the gradient < sfgrmin
[1] "Running with  175  repeats"
Square root transformation

rarified_count <- rrarefy(t(count_input),30)

#Convert the input into physeq objects
rare_count_phy <- otu_table(t(rarified_count), taxa_are_rows=T)
sample_info_tab_phy <- sample_data(duplicated_info)
rare_physeq <- phyloseq(rare_count_phy,sample_info_tab_phy)
  
#Perform the Ordination calculation
vst_pcoa <- ordinate(rare_physeq, method="NMDS", distance="bray")

Wisconsin double standardization
Run 0 stress 0.203313 
Run 1 stress 0.1998393 
... New best solution
... Procrustes: rmse 0.06275844  max resid 0.2040227 
Run 2 stress 0.1974667 
... New best solution
... Procrustes: rmse 0.06213639  max resid 0.254214 
Run 3 stress 0.2083041 
Run 4 stress 0.1974495 
... New best solution
... Procrustes: rmse 0.1244239  max resid 0.2606505 
Run 5 stress 0.1956229 
... New best solution
... Procrustes: rmse 0.1219681  max resid 0.2345158 
Run 6 stress 0.1987877 
Run 7 stress 0.1974475 
Run 8 stress 0.1986037 
Run 9 stress 0.2033898 
Run 10 stress 0.1962578 
Run 11 stress 0.1961679 
Run 12 stress 0.2014072 
Run 13 stress 0.1996942 
Run 14 stress 0.1999438 
Run 15 stress 0.2041151 
Run 16 stress 0.1951754 
... New best solution
... Procrustes: rmse 0.01070483  max resid 0.04630114 
Run 17 stress 0.2006656 
Run 18 stress 0.2082213 
Run 19 stress 0.1996263 
Run 20 stress 0.2006656 
*** Best solution was not repeated -- monoMDS stopping criteria:
     1: no. of iterations >= maxit
    18: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

  
#Convert ordination result into a data fram object
my_plot2 <- (plot_ordination(rare_physeq, vst_pcoa,justDF = TRUE))

just_positions <- my_plot2[,1:2]

cluster_output <- kmeans(just_positions,3, iter.max = 10, nstart = 25)


cluster_est <- cluster_output$cluster
cluster_true_pre <- sample_info_tab$location

cluster_names <- list()
cluster_true <- cluster_true_pre
for (x in 1:length(cluster_true_pre)) {
  if (!(cluster_true_pre[x] %in% cluster_names)){
    cluster_names <- append(cluster_names,cluster_true_pre[x])
  }
  cluster_true[x] <- as.integer(which(cluster_names == cluster_true_pre[x]))
  cluster_true <- as.numeric(unlist(cluster_true))
  }

Warning: NAs introduced by coercion

rand.index(cluster_est, cluster_true)

[1] 0.745098

plot_ordination(rare_physeq, vst_pcoa, color="location", shape="type")

cluster_est <- cluster_output$cluster
cluster_true_pre <- sample_info_tab$location

cluster_names <- list()
cluster_true <- cluster_true_pre
for (x in 1:length(cluster_true_pre)) {
  if (!(cluster_true_pre[x] %in% cluster_names)){
    cluster_names <- append(cluster_names,cluster_true_pre[x])
  }
  cluster_true[x] <- as.integer(which(cluster_names == cluster_true_pre[x]))
  cluster_true <- as.numeric(unlist(cluster_true))
  }

Warning: NAs introduced by coercion

rand.index(cluster_est, cluster_true)

[1] 0.7860963

Rand Index calculation (extrinsic similarity)

complete_ordination <- repeat_raref(count_data, sample_info_tab, repeats = 10, threshold = 1000, "sample_id", "location")

Square root transformation
Wisconsin double standardization
Run 0 stress 0.2159642 
Run 1 stress 0.2150414 
... New best solution
... Procrustes: rmse 0.01149688  max resid 0.09973348 
Run 2 stress 0.2280966 
Run 3 stress 0.2078487 
... New best solution
... Procrustes: rmse 0.01582487  max resid 0.1029343 
Run 4 stress 0.2688891 
Run 5 stress 0.2054215 
... New best solution
... Procrustes: rmse 0.007557131  max resid 0.0635759 
Run 6 stress 0.2098344 
Run 7 stress 0.2054215 
... Procrustes: rmse 9.725294e-06  max resid 5.587627e-05 
... Similar to previous best
Run 8 stress 0.2134476 
Run 9 stress 0.2462353 
Run 10 stress 0.2080794 
Run 11 stress 0.2071042 
Run 12 stress 0.2478765 
Run 13 stress 0.2432921 
Run 14 stress 0.2298165 
Run 15 stress 0.234119 
Run 16 stress 0.206836 
Run 17 stress 0.2628698 
Run 18 stress 0.2080089 
Run 19 stress 0.2528086 
Run 20 stress 0.2321972 
*** Best solution repeated 1 times

just_positions <- complete_ordination[,1:2]
cluster_output <- kmeans(just_positions,3, iter.max = 10, nstart = 25)


cluster_est <- cluster_output$cluster
cluster_true_pre <- complete_ordination$location

cluster_names <- list()
cluster_true <- cluster_true_pre
for (x in 1:length(cluster_true_pre)) {
  if (!(cluster_true_pre[x] %in% cluster_names)){
    cluster_names <- append(cluster_names,cluster_true_pre[x])
  }
  cluster_true[x] <- as.integer(which(cluster_names == cluster_true_pre[x]))
  cluster_true <- as.numeric(unlist(cluster_true))
  }

Warning: NAs introduced by coercion

rand.index(cluster_est, cluster_true)

[1] 0.8300191

FM_index(cluster_est, cluster_true)

[1] 0.7537729
attr(,"E_FM")
[1] 0.3442455
attr(,"V_FM")
[1] 7.157823e-06

#library(clusterSim)

index_value <- function(repeats, threshold) {
  complete_ordination <- repeat_raref(count_data, sample_info_tab, repeats = repeats, threshold = threshold, "sample_id", "location")
  
  
  just_positions <- complete_ordination[,1:2]
  
  
  #Create matrix with mean position for each sample
  sum_test <- as.data.frame(rowsum(just_positions, complete_ordination$sample_id))
  scale_factor <- length(t(just_positions))/length(t(sum_test))
  sum_test <- as.data.frame(sum_test/scale_factor)
  
  
  cluster_true_pre <- complete_ordination$location[1:(length(t(sum_test))/2)]
  
  cluster_names <- list()
  cluster_true <- cluster_true_pre
  # for (x in 1:length(cluster_true_pre)) {
  #   if (!(cluster_true_pre[x] %in% cluster_names)){
  #     cluster_names <- append(cluster_names,cluster_true_pre[x])
  #   }
  #   cluster_true[x] <- as.integer(which(cluster_names == cluster_true_pre[x]))
  #   cluster_true <- as.numeric(unlist(cluster_true))
  #   }
  
  for (x in 1:length(cluster_true_pre)) {
    if (cluster_true_pre[x] == "VK5") {
      cluster_true[x] <- 1
    }else {
      cluster_true[x] <-2
    }
  }
  
  cluster_true <- as.numeric(unlist(cluster_true))
  return(index.G1(sum_test,cluster_true))
  
}

index_value_full <- function(repeats, threshold) {
  complete_ordination <- repeat_raref(count_data, sample_info_tab, repeats = repeats, threshold = threshold, "NMDS", "location", shapeb = "type")
  
  
  just_positions <- complete_ordination[,1:2]
  cluster_true_pre <- complete_ordination$location
  
  cluster_names <- list()
  cluster_true <- cluster_true_pre
  # for (x in 1:length(cluster_true_pre)) {
  #   if (!(cluster_true_pre[x] %in% cluster_names)){
  #     cluster_names <- append(cluster_names,cluster_true_pre[x])
  #   }
  #   cluster_true[x] <- as.integer(which(cluster_names == cluster_true_pre[x]))
  #   cluster_true <- as.numeric(unlist(cluster_true))
  #   }
  
  for (x in 1:length(cluster_true_pre)) {
    if (cluster_true_pre[x] == "VK5") {
      cluster_true[x] <- 1
    }else {
      cluster_true[x] <-2
    }
  }
  
  cluster_true <- as.numeric(unlist(cluster_true))
  return(index.G1(just_positions,cluster_true))
  
}

#repeat_list <- c(1, 10, 25, 50)
threshhold <- c(10, 10, 10, 20, 20, 20, 30,30,30, 40,40,40, 50,50,50, 75,75,75, 100,100,100, 150,150,150, 200,200,200)
#index_data <- matrix(ncol = 3, nrow = 0)
colnames(index_data) <- c("Repeat Amount", "Threshold", "Index")


for (x in threshhold) {
  print(paste("Running with ",x," threshold"))
  index <- index_value_full(20, x)
  index_data <- rbind(index_data, c("20", x, index))
}

[1] "Running with  10  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1879065 
Run 1 stress 0.2351654 
Run 2 stress 0.2279338 
Run 3 stress 0.1947221 
Run 4 stress 0.2117983 
Run 5 stress 0.1926417 
Run 6 stress 0.187908 
... Procrustes: rmse 0.0009576278  max resid 0.009528393 
... Similar to previous best
Run 7 stress 0.1891428 
Run 8 stress 0.2017835 
Run 9 stress 0.2157631 
Run 10 stress 0.1899759 
Run 11 stress 0.1899244 
Run 12 stress 0.1888526 
Run 13 stress 0.1888713 
Run 14 stress 0.1931819 
Run 15 stress 0.1946292 
Run 16 stress 0.196437 
Run 17 stress 0.1913327 
Run 18 stress 0.2105558 
Run 19 stress 0.1960642 
Run 20 stress 0.1884513 
*** Best solution repeated 1 times
[1] "Running with  10  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1901362 
Run 1 stress 0.1901082 
... New best solution
... Procrustes: rmse 0.01463621  max resid 0.2203599 
Run 2 stress 0.2042342 
Run 3 stress 0.1970064 
Run 4 stress 0.1959539 
Run 5 stress 0.1939417 
Run 6 stress 0.2116451 
Run 7 stress 0.1918994 
Run 8 stress 0.1995829 
Run 9 stress 0.2124489 
Run 10 stress 0.2097084 
Run 11 stress 0.2040889 
Run 12 stress 0.1904782 
... Procrustes: rmse 0.01167574  max resid 0.211441 
Run 13 stress 0.1926301 
Run 14 stress 0.1986335 
Run 15 stress 0.1992735 
Run 16 stress 0.198989 
Run 17 stress 0.1930341 
Run 18 stress 0.19608 
Run 19 stress 0.1933617 
Run 20 stress 0.1988223 
*** Best solution was not repeated -- monoMDS stopping criteria:
    19: no. of iterations >= maxit
     1: scale factor of the gradient < sfgrmin
[1] "Running with  10  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1912192 
Run 1 stress 0.2080146 
Run 2 stress 0.2225626 
Run 3 stress 0.1943057 
Run 4 stress 0.2174577 
Run 5 stress 0.2004843 
Run 6 stress 0.1928852 
Run 7 stress 0.1926066 
Run 8 stress 0.1902666 
... New best solution
... Procrustes: rmse 0.01332018  max resid 0.2244647 
Run 9 stress 0.1918422 
Run 10 stress 0.2152146 
Run 11 stress 0.2135304 
Run 12 stress 0.1983539 
Run 13 stress 0.196493 
Run 14 stress 0.2046752 
Run 15 stress 0.1930336 
Run 16 stress 0.2356528 
Run 17 stress 0.1947518 
Run 18 stress 0.1936484 
Run 19 stress 0.1971467 
Run 20 stress 0.194448 
*** Best solution was not repeated -- monoMDS stopping criteria:
    18: no. of iterations >= maxit
     2: stress ratio > sratmax
[1] "Running with  20  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.208276 
Run 1 stress 0.2010141 
... New best solution
... Procrustes: rmse 0.01394661  max resid 0.16227 
Run 2 stress 0.2075086 
Run 3 stress 0.2468217 
Run 4 stress 0.2040368 
Run 5 stress 0.2255259 
Run 6 stress 0.2009646 
... New best solution
... Procrustes: rmse 0.01093212  max resid 0.1472212 
Run 7 stress 0.2074124 
Run 8 stress 0.2070586 
Run 9 stress 0.2044067 
Run 10 stress 0.2078921 
Run 11 stress 0.2062613 
Run 12 stress 0.2113081 
Run 13 stress 0.2323731 
Run 14 stress 0.1996351 
... New best solution
... Procrustes: rmse 0.008592032  max resid 0.1476538 
Run 15 stress 0.2227702 
Run 16 stress 0.2215638 
Run 17 stress 0.2307306 
Run 18 stress 0.2228376 
Run 19 stress 0.2015223 
Run 20 stress 0.2064459 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
    14: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin
[1] "Running with  20  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2092519 
Run 1 stress 0.2052065 
... New best solution
... Procrustes: rmse 0.01733205  max resid 0.1878514 
Run 2 stress 0.2097542 
Run 3 stress 0.2326529 
Run 4 stress 0.2445815 
Run 5 stress 0.203301 
... New best solution
... Procrustes: rmse 0.01330983  max resid 0.1942919 
Run 6 stress 0.2046527 
Run 7 stress 0.2008142 
... New best solution
... Procrustes: rmse 0.01096243  max resid 0.1969584 
Run 8 stress 0.22198 
Run 9 stress 0.2053797 
Run 10 stress 0.2008129 
... New best solution
... Procrustes: rmse 0.00012504  max resid 0.00212188 
... Similar to previous best
Run 11 stress 0.2191037 
Run 12 stress 0.2525629 
Run 13 stress 0.2054609 
Run 14 stress 0.2026013 
Run 15 stress 0.2082459 
Run 16 stress 0.2029501 
Run 17 stress 0.2219549 
Run 18 stress 0.2300808 
Run 19 stress 0.2024969 
Run 20 stress 0.2030065 
*** Best solution repeated 1 times
[1] "Running with  20  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1974386 
Run 1 stress 0.2234282 
Run 2 stress 0.2311522 
Run 3 stress 0.202053 
Run 4 stress 0.2211812 
Run 5 stress 0.2125345 
Run 6 stress 0.2007728 
Run 7 stress 0.2064168 
Run 8 stress 0.2135481 
Run 9 stress 0.2026546 
Run 10 stress 0.2006983 
Run 11 stress 0.2009266 
Run 12 stress 0.2225039 
Run 13 stress 0.1992542 
Run 14 stress 0.1974211 
... New best solution
... Procrustes: rmse 0.01092241  max resid 0.2092013 
Run 15 stress 0.1974231 
... Procrustes: rmse 0.0005328633  max resid 0.008370599 
... Similar to previous best
Run 16 stress 0.2004732 
Run 17 stress 0.1984804 
Run 18 stress 0.1989647 
Run 19 stress 0.1998359 
Run 20 stress 0.2225625 
*** Best solution repeated 1 times
[1] "Running with  30  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2121851 
Run 1 stress 0.2055815 
... New best solution
... Procrustes: rmse 0.01595801  max resid 0.1600952 
Run 2 stress 0.2308756 
Run 3 stress 0.2043292 
... New best solution
... Procrustes: rmse 0.01330054  max resid 0.1752411 
Run 4 stress 0.2298599 
Run 5 stress 0.233246 
Run 6 stress 0.2136824 
Run 7 stress 0.2130154 
Run 8 stress 0.2071256 
Run 9 stress 0.2102505 
Run 10 stress 0.2294925 
Run 11 stress 0.2103372 
Run 12 stress 0.2044111 
... Procrustes: rmse 0.01359241  max resid 0.175295 
Run 13 stress 0.2316389 
Run 14 stress 0.2440312 
Run 15 stress 0.2110576 
Run 16 stress 0.2043978 
... Procrustes: rmse 0.01194033  max resid 0.1777093 
Run 17 stress 0.2319735 
Run 18 stress 0.2031971 
... New best solution
... Procrustes: rmse 0.01089643  max resid 0.1775277 
Run 19 stress 0.2051502 
Run 20 stress 0.2553234 
*** Best solution was not repeated -- monoMDS stopping criteria:
     7: no. of iterations >= maxit
    11: stress ratio > sratmax
     2: scale factor of the gradient < sfgrmin
[1] "Running with  30  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2032955 
Run 1 stress 0.2266979 
Run 2 stress 0.2310152 
Run 3 stress 0.2051493 
Run 4 stress 0.231805 
Run 5 stress 0.2055934 
Run 6 stress 0.231565 
Run 7 stress 0.2338786 
Run 8 stress 0.213489 
Run 9 stress 0.2351557 
Run 10 stress 0.2039548 
Run 11 stress 0.2147903 
Run 12 stress 0.236825 
Run 13 stress 0.2207542 
Run 14 stress 0.2299267 
Run 15 stress 0.2050121 
Run 16 stress 0.2038053 
Run 17 stress 0.2237606 
Run 18 stress 0.2298634 
Run 19 stress 0.23732 
Run 20 stress 0.2510337 
*** Best solution was not repeated -- monoMDS stopping criteria:
     8: no. of iterations >= maxit
     8: stress ratio > sratmax
     4: scale factor of the gradient < sfgrmin
[1] "Running with  30  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2022775 
Run 1 stress 0.2069794 
Run 2 stress 0.2100063 
Run 3 stress 0.2271791 
Run 4 stress 0.2114856 
Run 5 stress 0.2511749 
Run 6 stress 0.2421795 
Run 7 stress 0.229245 
Run 8 stress 0.2036348 
Run 9 stress 0.2253002 
Run 10 stress 0.2076309 
Run 11 stress 0.2374795 
Run 12 stress 0.211179 
Run 13 stress 0.2083897 
Run 14 stress 0.2049113 
Run 15 stress 0.2300812 
Run 16 stress 0.2274549 
Run 17 stress 0.2054203 
Run 18 stress 0.239142 
Run 19 stress 0.2207944 
Run 20 stress 0.2312561 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
    11: stress ratio > sratmax
     4: scale factor of the gradient < sfgrmin
[1] "Running with  40  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2103896 
Run 1 stress 0.2190143 
Run 2 stress 0.2366566 
Run 3 stress 0.2370157 
Run 4 stress 0.2559342 
Run 5 stress 0.209979 
... New best solution
... Procrustes: rmse 0.002232265  max resid 0.01167297 
Run 6 stress 0.2335863 
Run 7 stress 0.2474291 
Run 8 stress 0.2264526 
Run 9 stress 0.2101734 
... Procrustes: rmse 0.007913271  max resid 0.1591419 
Run 10 stress 0.2414867 
Run 11 stress 0.2121263 
Run 12 stress 0.2378325 
Run 13 stress 0.2361034 
Run 14 stress 0.2190178 
Run 15 stress 0.2131337 
Run 16 stress 0.2379373 
Run 17 stress 0.2391365 
Run 18 stress 0.2303841 
Run 19 stress 0.2440825 
Run 20 stress 0.2346862 
*** Best solution was not repeated -- monoMDS stopping criteria:
     9: no. of iterations >= maxit
    10: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin
[1] "Running with  40  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2148464 
Run 1 stress 0.2472881 
Run 2 stress 0.2087489 
... New best solution
... Procrustes: rmse 0.01591775  max resid 0.1668105 
Run 3 stress 0.2273732 
Run 4 stress 0.2087817 
... Procrustes: rmse 0.01490556  max resid 0.1684122 
Run 5 stress 0.232904 
Run 6 stress 0.2396605 
Run 7 stress 0.2464808 
Run 8 stress 0.2309457 
Run 9 stress 0.2093931 
Run 10 stress 0.2329445 
Run 11 stress 0.2142143 
Run 12 stress 0.2342316 
Run 13 stress 0.2088601 
... Procrustes: rmse 0.01470783  max resid 0.1674293 
Run 14 stress 0.2085718 
... New best solution
... Procrustes: rmse 0.0129591  max resid 0.1698823 
Run 15 stress 0.235456 
Run 16 stress 0.2365738 
Run 17 stress 0.2068191 
... New best solution
... Procrustes: rmse 0.01343894  max resid 0.1705109 
Run 18 stress 0.2558904 
Run 19 stress 0.2452189 
Run 20 stress 0.2093018 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
     9: stress ratio > sratmax
     9: scale factor of the gradient < sfgrmin
[1] "Running with  40  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.207224 
Run 1 stress 0.2323069 
Run 2 stress 0.2346807 
Run 3 stress 0.2327263 
Run 4 stress 0.2096719 
Run 5 stress 0.2311162 
Run 6 stress 0.2567539 
Run 7 stress 0.2327764 
Run 8 stress 0.216323 
Run 9 stress 0.2380783 
Run 10 stress 0.2409997 
Run 11 stress 0.2250466 
Run 12 stress 0.253229 
Run 13 stress 0.2153763 
Run 14 stress 0.2292307 
Run 15 stress 0.2337005 
Run 16 stress 0.2318753 
Run 17 stress 0.2302579 
Run 18 stress 0.2169229 
Run 19 stress 0.2097838 
Run 20 stress 0.2317223 
*** Best solution was not repeated -- monoMDS stopping criteria:
     4: no. of iterations >= maxit
     9: stress ratio > sratmax
     7: scale factor of the gradient < sfgrmin
[1] "Running with  50  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2139669 
Run 1 stress 0.2100292 
... New best solution
... Procrustes: rmse 0.01670073  max resid 0.1521534 
Run 2 stress 0.2315663 
Run 3 stress 0.2257774 
Run 4 stress 0.2093604 
... New best solution
... Procrustes: rmse 0.01025766  max resid 0.1585246 
Run 5 stress 0.2394295 
Run 6 stress 0.240939 
Run 7 stress 0.2322853 
Run 8 stress 0.2250766 
Run 9 stress 0.2442281 
Run 10 stress 0.228365 
Run 11 stress 0.2474384 
Run 12 stress 0.2057771 
... New best solution
... Procrustes: rmse 0.009738888  max resid 0.1453372 
Run 13 stress 0.2335074 
Run 14 stress 0.2070094 
Run 15 stress 0.2335297 
Run 16 stress 0.2251606 
Run 17 stress 0.214186 
Run 18 stress 0.2329643 
Run 19 stress 0.2271108 
Run 20 stress 0.2242644 
*** Best solution was not repeated -- monoMDS stopping criteria:
     6: no. of iterations >= maxit
    10: stress ratio > sratmax
     4: scale factor of the gradient < sfgrmin
[1] "Running with  50  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2156402 
Run 1 stress 0.2300636 
Run 2 stress 0.2284606 
Run 3 stress 0.2100246 
... New best solution
... Procrustes: rmse 0.01832015  max resid 0.1790795 
Run 4 stress 0.2187009 
Run 5 stress 0.226346 
Run 6 stress 0.2533131 
Run 7 stress 0.2117277 
Run 8 stress 0.2347731 
Run 9 stress 0.2069812 
... New best solution
... Procrustes: rmse 0.01169  max resid 0.1861915 
Run 10 stress 0.2184447 
Run 11 stress 0.2147922 
Run 12 stress 0.2065306 
... New best solution
... Procrustes: rmse 0.007204208  max resid 0.1634066 
Run 13 stress 0.2101619 
Run 14 stress 0.2105052 
Run 15 stress 0.2373002 
Run 16 stress 0.20695 
... Procrustes: rmse 0.007190924  max resid 0.1700698 
Run 17 stress 0.2448715 
Run 18 stress 0.2069518 
... Procrustes: rmse 0.007383915  max resid 0.1677263 
Run 19 stress 0.2465397 
Run 20 stress 0.2344034 
*** Best solution was not repeated -- monoMDS stopping criteria:
     6: no. of iterations >= maxit
     9: stress ratio > sratmax
     5: scale factor of the gradient < sfgrmin
[1] "Running with  50  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2146762 
Run 1 stress 0.2393342 
Run 2 stress 0.2394557 
Run 3 stress 0.2157252 
Run 4 stress 0.2152456 
Run 5 stress 0.2515705 
Run 6 stress 0.2434822 
Run 7 stress 0.2484161 
Run 8 stress 0.2104424 
... New best solution
... Procrustes: rmse 0.01256799  max resid 0.1548368 
Run 9 stress 0.2368677 
Run 10 stress 0.2386055 
Run 11 stress 0.2433088 
Run 12 stress 0.217088 
Run 13 stress 0.2440282 
Run 14 stress 0.2328256 
Run 15 stress 0.2088121 
... New best solution
... Procrustes: rmse 0.01118308  max resid 0.1559095 
Run 16 stress 0.2328368 
Run 17 stress 0.2142142 
Run 18 stress 0.2139738 
Run 19 stress 0.2109498 
Run 20 stress 0.2099373 
*** Best solution was not repeated -- monoMDS stopping criteria:
     1: no. of iterations >= maxit
     8: stress ratio > sratmax
    11: scale factor of the gradient < sfgrmin
[1] "Running with  75  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2224438 
Run 1 stress 0.2166283 
... New best solution
... Procrustes: rmse 0.01607284  max resid 0.1347807 
Run 2 stress 0.2375431 
Run 3 stress 0.2411074 
Run 4 stress 0.2362093 
Run 5 stress 0.2294638 
Run 6 stress 0.2114948 
... New best solution
... Procrustes: rmse 0.01068947  max resid 0.1395129 
Run 7 stress 0.2339678 
Run 8 stress 0.2629333 
Run 9 stress 0.2280436 
Run 10 stress 0.2399913 
Run 11 stress 0.2369366 
Run 12 stress 0.2487186 
Run 13 stress 0.2332686 
Run 14 stress 0.2395857 
Run 15 stress 0.2174297 
Run 16 stress 0.227616 
Run 17 stress 0.245257 
Run 18 stress 0.2537805 
Run 19 stress 0.2641542 
Run 20 stress 0.2128824 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
     6: stress ratio > sratmax
     9: scale factor of the gradient < sfgrmin
[1] "Running with  75  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.212657 
Run 1 stress 0.2163426 
Run 2 stress 0.23807 
Run 3 stress 0.222213 
Run 4 stress 0.2141724 
Run 5 stress 0.2524655 
Run 6 stress 0.2528646 
Run 7 stress 0.2501903 
Run 8 stress 0.2404367 
Run 9 stress 0.2428561 
Run 10 stress 0.2129548 
... Procrustes: rmse 0.01280873  max resid 0.1480826 
Run 11 stress 0.2115998 
... New best solution
... Procrustes: rmse 0.01041513  max resid 0.1496079 
Run 12 stress 0.2110805 
... New best solution
... Procrustes: rmse 0.008034991  max resid 0.1512215 
Run 13 stress 0.2189774 
Run 14 stress 0.2317403 
Run 15 stress 0.2132136 
Run 16 stress 0.2428716 
Run 17 stress 0.2393786 
Run 18 stress 0.2502293 
Run 19 stress 0.2331742 
Run 20 stress 0.2474184 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
     6: stress ratio > sratmax
    12: scale factor of the gradient < sfgrmin
[1] "Running with  75  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2129167 
Run 1 stress 0.2367609 
Run 2 stress 0.2202147 
Run 3 stress 0.235514 
Run 4 stress 0.2133565 
... Procrustes: rmse 0.01026145  max resid 0.1295895 
Run 5 stress 0.2242016 
Run 6 stress 0.2143594 
Run 7 stress 0.2221496 
Run 8 stress 0.2202375 
Run 9 stress 0.2122574 
... New best solution
... Procrustes: rmse 0.01233678  max resid 0.1286866 
Run 10 stress 0.2396891 
Run 11 stress 0.2223138 
Run 12 stress 0.245494 
Run 13 stress 0.2155583 
Run 14 stress 0.2374963 
Run 15 stress 0.2621461 
Run 16 stress 0.2218124 
Run 17 stress 0.2167978 
Run 18 stress 0.2605603 
Run 19 stress 0.2194657 
Run 20 stress 0.2323379 
*** Best solution was not repeated -- monoMDS stopping criteria:
     6: no. of iterations >= maxit
     5: stress ratio > sratmax
     9: scale factor of the gradient < sfgrmin
[1] "Running with  100  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2221003 
Run 1 stress 0.2417038 
Run 2 stress 0.2515101 
Run 3 stress 0.2280799 
Run 4 stress 0.2246659 
Run 5 stress 0.2340205 
Run 6 stress 0.237241 
Run 7 stress 0.2183764 
... New best solution
... Procrustes: rmse 0.0151979  max resid 0.1251456 
Run 8 stress 0.2176732 
... New best solution
... Procrustes: rmse 0.01269477  max resid 0.1377417 
Run 9 stress 0.2268045 
Run 10 stress 0.2473608 
Run 11 stress 0.2134028 
... New best solution
... Procrustes: rmse 0.0114971  max resid 0.1389168 
Run 12 stress 0.2399769 
Run 13 stress 0.2395119 
Run 14 stress 0.2515725 
Run 15 stress 0.214438 
Run 16 stress 0.2339173 
Run 17 stress 0.2361692 
Run 18 stress 0.2505785 
Run 19 stress 0.2635414 
Run 20 stress 0.2241874 
*** Best solution was not repeated -- monoMDS stopping criteria:
     7: no. of iterations >= maxit
     6: stress ratio > sratmax
     7: scale factor of the gradient < sfgrmin
[1] "Running with  100  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2125345 
Run 1 stress 0.2403776 
Run 2 stress 0.2443304 
Run 3 stress 0.2394101 
Run 4 stress 0.2382654 
Run 5 stress 0.2397149 
Run 6 stress 0.2189959 
Run 7 stress 0.2402836 
Run 8 stress 0.2405579 
Run 9 stress 0.2382527 
Run 10 stress 0.214554 
Run 11 stress 0.2392091 
Run 12 stress 0.2300048 
Run 13 stress 0.2367911 
Run 14 stress 0.2179475 
Run 15 stress 0.2241107 
Run 16 stress 0.2199652 
Run 17 stress 0.2165537 
Run 18 stress 0.2213124 
Run 19 stress 0.2416903 
Run 20 stress 0.2537032 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
     7: stress ratio > sratmax
    11: scale factor of the gradient < sfgrmin
[1] "Running with  100  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2227097 
Run 1 stress 0.2473109 
Run 2 stress 0.2118124 
... New best solution
... Procrustes: rmse 0.01590452  max resid 0.1427848 
Run 3 stress 0.2415765 
Run 4 stress 0.2210779 
Run 5 stress 0.2454134 
Run 6 stress 0.2192524 
Run 7 stress 0.2123718 
Run 8 stress 0.2166801 
Run 9 stress 0.2136508 
Run 10 stress 0.2572876 
Run 11 stress 0.2128444 
Run 12 stress 0.2185548 
Run 13 stress 0.2537878 
Run 14 stress 0.2203629 
Run 15 stress 0.2408526 
Run 16 stress 0.2438652 
Run 17 stress 0.2215665 
Run 18 stress 0.2357806 
Run 19 stress 0.2363958 
Run 20 stress 0.2129073 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
     7: stress ratio > sratmax
     8: scale factor of the gradient < sfgrmin
[1] "Running with  150  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2234699 
Run 1 stress 0.2140673 
... New best solution
... Procrustes: rmse 0.01466064  max resid 0.1205849 
Run 2 stress 0.2336796 
Run 3 stress 0.2493266 
Run 4 stress 0.2455777 
Run 5 stress 0.2310956 
Run 6 stress 0.2407018 
Run 7 stress 0.2408338 
Run 8 stress 0.2383997 
Run 9 stress 0.2398284 
Run 10 stress 0.2495567 
Run 11 stress 0.2440914 
Run 12 stress 0.2627829 
Run 13 stress 0.2194716 
Run 14 stress 0.2352672 
Run 15 stress 0.2527393 
Run 16 stress 0.235836 
Run 17 stress 0.2486583 
Run 18 stress 0.2343089 
Run 19 stress 0.2485671 
Run 20 stress 0.2553643 
*** Best solution was not repeated -- monoMDS stopping criteria:
     4: no. of iterations >= maxit
     3: stress ratio > sratmax
    13: scale factor of the gradient < sfgrmin
[1] "Running with  150  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2243596 
Run 1 stress 0.226119 
Run 2 stress 0.2145907 
... New best solution
... Procrustes: rmse 0.0161654  max resid 0.1217609 
Run 3 stress 0.2487656 
Run 4 stress 0.2452631 
Run 5 stress 0.2536605 
Run 6 stress 0.2385913 
Run 7 stress 0.2386343 
Run 8 stress 0.2219692 
Run 9 stress 0.2197042 
Run 10 stress 0.2476442 
Run 11 stress 0.2375095 
Run 12 stress 0.2212111 
Run 13 stress 0.2356874 
Run 14 stress 0.2373525 
Run 15 stress 0.2150191 
... Procrustes: rmse 0.008840583  max resid 0.1165975 
Run 16 stress 0.2389051 
Run 17 stress 0.2233348 
Run 18 stress 0.2213175 
Run 19 stress 0.2446109 
Run 20 stress 0.2655289 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
     4: stress ratio > sratmax
    13: scale factor of the gradient < sfgrmin
[1] "Running with  150  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2240634 
Run 1 stress 0.2241869 
... Procrustes: rmse 0.01982327  max resid 0.1300086 
Run 2 stress 0.2408061 
Run 3 stress 0.2213032 
... New best solution
... Procrustes: rmse 0.01941456  max resid 0.1305445 
Run 4 stress 0.2530636 
Run 5 stress 0.2404874 
Run 6 stress 0.2434094 
Run 7 stress 0.2220371 
Run 8 stress 0.2498916 
Run 9 stress 0.2346629 
Run 10 stress 0.2505836 
Run 11 stress 0.2407792 
Run 12 stress 0.2484317 
Run 13 stress 0.2559916 
Run 14 stress 0.2452553 
Run 15 stress 0.2147341 
... New best solution
... Procrustes: rmse 0.01119022  max resid 0.1342018 
Run 16 stress 0.2447422 
Run 17 stress 0.238059 
Run 18 stress 0.2411647 
Run 19 stress 0.2435643 
Run 20 stress 0.2396175 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
     6: stress ratio > sratmax
    12: scale factor of the gradient < sfgrmin
[1] "Running with  200  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.224605 
Run 1 stress 0.2664843 
Run 2 stress 0.2183212 
... New best solution
... Procrustes: rmse 0.01452049  max resid 0.1122702 
Run 3 stress 0.2514294 
Run 4 stress 0.2129609 
... New best solution
... Procrustes: rmse 0.008589928  max resid 0.1215177 
Run 5 stress 0.216691 
Run 6 stress 0.222422 
Run 7 stress 0.2367106 
Run 8 stress 0.2383715 
Run 9 stress 0.2376185 
Run 10 stress 0.2544464 
Run 11 stress 0.2440541 
Run 12 stress 0.2144243 
Run 13 stress 0.2394312 
Run 14 stress 0.2665107 
Run 15 stress 0.2241192 
Run 16 stress 0.2427227 
Run 17 stress 0.2244602 
Run 18 stress 0.2230144 
Run 19 stress 0.2403596 
Run 20 stress 0.216699 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
     1: stress ratio > sratmax
    17: scale factor of the gradient < sfgrmin
[1] "Running with  200  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.223739 
Run 1 stress 0.2397261 
Run 2 stress 0.2118603 
... New best solution
... Procrustes: rmse 0.01602588  max resid 0.1237196 
Run 3 stress 0.241145 
Run 4 stress 0.2364524 
Run 5 stress 0.212158 
... Procrustes: rmse 0.008621638  max resid 0.1288537 
Run 6 stress 0.2405229 
Run 7 stress 0.2321443 
Run 8 stress 0.2394778 
Run 9 stress 0.2516443 
Run 10 stress 0.2148492 
Run 11 stress 0.2359846 
Run 12 stress 0.4199394 
Run 13 stress 0.2408343 
Run 14 stress 0.2598164 
Run 15 stress 0.2413239 
Run 16 stress 0.2372033 
Run 17 stress 0.2139371 
Run 18 stress 0.2393511 
Run 19 stress 0.2153808 
Run 20 stress 0.2294073 
*** Best solution was not repeated -- monoMDS stopping criteria:
     6: stress ratio > sratmax
    14: scale factor of the gradient < sfgrmin
[1] "Running with  200  threshold"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.222938 
Run 1 stress 0.2385912 
Run 2 stress 0.2393115 
Run 3 stress 0.2125013 
... New best solution
... Procrustes: rmse 0.01686376  max resid 0.1368359 
Run 4 stress 0.2141032 
Run 5 stress 0.2506911 
Run 6 stress 0.2126779 
... Procrustes: rmse 0.009411201  max resid 0.1419478 
Run 7 stress 0.2135901 
Run 8 stress 0.2488331 
Run 9 stress 0.2384335 
Run 10 stress 0.24222 
Run 11 stress 0.2647757 
Run 12 stress 0.2419426 
Run 13 stress 0.2296467 
Run 14 stress 0.2245041 
Run 15 stress 0.2604258 
Run 16 stress 0.2529661 
Run 17 stress 0.2403464 
Run 18 stress 0.244037 
Run 19 stress 0.2438353 
Run 20 stress 0.2404865 
*** Best solution was not repeated -- monoMDS stopping criteria:
     4: no. of iterations >= maxit
     3: stress ratio > sratmax
    13: scale factor of the gradient < sfgrmin

plot(x = index_data[,2], y = index_data[,3])

colnames(index_data) <- c("Repeat_Amount", "Threshold", "Index")
index_data <- as.data.frame(index_data)
index_data$Threshold <- as.numeric(index_data$Threshold)
index_data$Index <- as.numeric(index_data$Index)
 
ggplot(index_data,aes(x=Threshold,y=Index,color=Repeat_Amount)) + geom_point() + labs(x = "Rarefaction Threshold", y ="Calinski-Harabasz pseudo F-statistic")

index_data_wor <- index_data[index_data$Repeat_Amount == "1",,]
ggplot(index_data_wor,aes(x=Threshold,y=Index,color=Repeat_Amount)) + geom_point() + labs(x = "Rarefaction Threshold", y ="Calinski-Harabasz pseudo F-statistic")

index <- index_value_full(20, 5)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1680435 
Run 1 stress 0.166982 
... New best solution
... Procrustes: rmse 0.02041138  max resid 0.2835901 
Run 2 stress 0.1919353 
Run 3 stress 0.1679044 
Run 4 stress 0.1685771 
Run 5 stress 0.1671102 
... Procrustes: rmse 0.01767279  max resid 0.2412741 
Run 6 stress 0.1711792 
Run 7 stress 0.190869 
Run 8 stress 0.1665096 
... New best solution
... Procrustes: rmse 0.01643726  max resid 0.2708281 
Run 9 stress 0.1683286 
Run 10 stress 0.1672181 
Run 11 stress 0.1664321 
... New best solution
... Procrustes: rmse 0.009156992  max resid 0.1821273 
Run 12 stress 0.1670267 
Run 13 stress 0.1672244 
Run 14 stress 0.167243 
Run 15 stress 0.1824932 
Run 16 stress 0.1680657 
Run 17 stress 0.169201 
Run 18 stress 0.1683534 
Run 19 stress 0.1669714 
Run 20 stress 0.1705786 
*** Best solution was not repeated -- monoMDS stopping criteria:
    20: no. of iterations >= maxit

index_data_r <- rbind(index_data_r, c("20", 5, index))

index_data_r$Index <- as.numeric(index_data_r$Index)
index_data_r$Threshold <- as.numeric(index_data_r$Threshold)
ggplot(index_data_r,aes(x=Threshold,y=Index,color=Repeat_Amount)) + geom_point() + labs(x = "Rarefaction Threshold", y ="Calinski-Harabasz pseudo F-statistic")